% This function solves for the the equivalent variation in the case of opt-out.
% Note: m1 is the amount by which income is augmented by the equalizing compensation (currently zero each time function is called), which can affect the contribution rate choice. In our final specifications, we keep the contribution rate fixed as we vary equalizing compensation (see footnote 24).

function[value]=ev_optout_int(m1,gamma,alpha,rho,t,m,maxmatch,maxcont,mm)

% Define "virtual income" above kink
Zbar=1+(1-t).*maxmatch.*(1+m).*(1-(1/(1+m)));

%Define EV in the case of opt out for various cases of x around the maximum
%match rate. 

%There are four cases, b/c x(theta,0) and x(theta,m1) can be above or below
%the match cap.

% Case 1: Both are below
% Case 2: Both are above
% Case 3: x*(theta,0) above and x*(theta,m1) below
% Case 4: x*(theta,0) below and x*(theta,m1) above

ev_ifxbelowM = exp(-gamma).*...
    (((xstar_ev(0,alpha,rho,t,m,Zbar,maxmatch,maxcont)+alpha)./...
    (xstar_ev(m1,alpha,rho,t,m,Zbar,maxmatch,maxcont)+alpha)).^rho).*...
    (1-((1-t)./(1+m)).*xstar_ev(0,alpha,rho,t,m,Zbar,maxmatch,maxcont))+(((1-t)./(1+m))...
    .*xstar_ev(m1,alpha,rho,t,m,Zbar,maxmatch,maxcont))- 1 ...
    +mm.*xstar_ev(m1,alpha,rho,t,m,Zbar,maxmatch,maxcont);

ev_ifxaboveM = exp(-gamma).*...
    (((xstar_ev(0,alpha,rho,t,m,Zbar,maxmatch,maxcont)+alpha)./...
    (xstar_ev(m1,alpha,rho,t,m,Zbar,maxmatch,maxcont)+alpha)).^rho).*...
    (Zbar-((1-t)).*xstar_ev(0,alpha,rho,t,m,Zbar,maxmatch,maxcont))+(((1-t))...
    .*xstar_ev(m1,alpha,rho,t,m,Zbar,maxmatch,maxcont))- Zbar ...
    +mm.*xstar_ev(m1,alpha,rho,t,m,Zbar,maxmatch,maxcont);

ev_case3 = exp(-gamma).*...
    (((xstar_ev(0,alpha,rho,t,m,Zbar,maxmatch,maxcont)+alpha)./...
    (xstar_ev(m1,alpha,rho,t,m,Zbar,maxmatch,maxcont)+alpha)).^rho).*...
    (Zbar-((1-t))*xstar_ev(0,alpha,rho,t,m,Zbar,maxmatch,maxcont))+((1-t)./(1+m))...
    .*xstar_ev(m1,alpha,rho,t,m,Zbar,maxmatch,maxcont)- 1 ...
    +mm.*xstar_ev(m1,alpha,rho,t,m,Zbar,maxmatch,maxcont);

ev_case4 = exp(-gamma).*...
    (((xstar_ev(0,alpha,rho,t,m,Zbar,maxmatch,maxcont)+alpha)./...
    (xstar_ev(m1,alpha,rho,t,m,Zbar,maxmatch,maxcont)+alpha)).^rho).*...
    (1-((1-t)./(1+m))*xstar_ev(0,alpha,rho,t,m,Zbar,maxmatch,maxcont))+((1-t)...
    .*xstar_ev(m1,alpha,rho,t,m,Zbar,maxmatch,maxcont))- Zbar ...
    +mm.*xstar_ev(m1,alpha,rho,t,m,Zbar,maxmatch,maxcont);

value = ev_ifxbelowM.*(xstar_ev(0,alpha,rho,t,m,Zbar,maxmatch,maxcont)<=maxmatch.*(1+m))...
    .*(xstar_ev(m1,alpha,rho,t,m,Zbar,maxmatch,maxcont)<=maxmatch.*(1+m))...
    + ev_ifxaboveM.*(xstar_ev(0,alpha,rho,t,m,Zbar,maxmatch,maxcont)>maxmatch.*(1+m))...
    .*(xstar_ev(m1,alpha,rho,t,m,Zbar,maxmatch,maxcont)>maxmatch*(1+m))...
    + ev_case3.*(xstar_ev(0,alpha,rho,t,m,Zbar,maxmatch,maxcont)>maxmatch.*(1+m))...
    .*(xstar_ev(m1,alpha,rho,t,m,Zbar,maxmatch,maxcont)<=maxmatch.*(1+m))...
    + ev_case4.*(xstar_ev(0,alpha,rho,t,m,Zbar,maxmatch,maxcont)<=maxmatch.*(1+m))...
    .*(xstar_ev(m1,alpha,rho,t,m,Zbar,maxmatch,maxcont)>maxmatch.*(1+m));

    % x*(theta,ev) function
    function[value]=xstar_ev(m1,alpha,rho,t,m,Zbar,maxmatch,maxcont)
    xstar_ifbelowM = (rho/(1+rho)).*((1+m).*(1+m1)./(1-t))-alpha./(1+rho);
    xstar_ifaboveM = (rho/(1+rho)).*((Zbar+m1)./(1-t))-alpha./(1+rho);
    xstar = xstar_ifbelowM.*(xstar_ifbelowM<=maxmatch.*(1+m)) + xstar_ifaboveM.*(xstar_ifaboveM>maxmatch.*(1+m))...
        + maxmatch.*(1+m).*(xstar_ifbelowM>maxmatch.*(1+m)).*(xstar_ifaboveM<=maxmatch.*(1+m));
    xstar=max(xstar,0);
    
    % Now, we account for the fact that the personal contribution,
    % x_personal, must be a discrete value (i.e. 1%, 2%, etc.)
    x_personal_opt=min(maxmatch.*(1+m),xstar)./(1+m)+max(0,xstar-(maxmatch.*(1+m)));
    
    % Calculate discrete % contributions above and below optimal personal
    % contribution. (Can't contribute < 0)
    x_personal_top=ceil(x_personal_opt.*100)./100;
    x_personal_top=max(x_personal_top,0);
    
    x_personal_bot=floor(x_personal_opt.*100)./100;
    x_personal_bot=max(x_personal_bot,0);
    
    % Can't contribute > maxcont
    x_personal_top=min(maxcont,x_personal_top);
    x_personal_bot=min(maxcont,x_personal_bot);
    
    % So the choice, x, would actually be one of these two discrete values
    x_top=x_personal_top+(m.*min(x_personal_top,maxmatch));
    x_bot=x_personal_bot+(m.*min(x_personal_bot,maxmatch));
    
    % The employee chooses the round number that gives the highest utility
    utility_top = ((rho.*log(x_top+alpha) + log(1+m1-((1-t)./(1+m)).*x_top)).*...
        (x_top<=maxmatch.*(1+m))... 
        +((rho.*log(x_top+alpha)+log(Zbar+m1-(1-t).*x_top)).*((x_top>maxmatch.*(1+m)))));
    utility_bot = ((rho.*log(x_bot+alpha) + log(1+m1-((1-t)./(1+m)).*x_bot)).*...
        (x_bot<=maxmatch.*(1+m)))... 
        +((rho.*log(x_bot+alpha)+log(Zbar+m1-(1-t).*x_bot)).*((x_bot>maxmatch.*(1+m))));
    
    value=(xstar>=0).*( (x_top).*(utility_top>=utility_bot) + (x_bot).*(utility_top<utility_bot));
    end

end